

clear


%% data
load('data_5y.mat')


%% pre-1950 data (excluding WWs)
d0=data((data.year>=1875&data.year<=1950)&isfinite(data.y),:);
d0(d0.year>=1914&d0.year<=1918,:)=[];
d0(d0.year>=1939&d0.year<=1945,:)=[];
ccodes=unique(d0.ccode);
countryFE=zeros(size(d0,1),length(ccodes));
for m=1:length(unique(ccodes))
    countryFE(:,m)=d0.ccode==ccodes(m);
end
years=unique(d0.year);
yearFE=zeros(size(d0,1),length(years));
for m=1:length(unique(years))
    yearFE(:,m)=d0.year==years(m);
end
avgGR=(d0.y'*countryFE)./sum(countryFE);


%% \theta
th0=0;
om_th=1;


%% \beta
b0A=mean(d0.y(d0.D==0));
b0D=mean(d0.y(d0.D==1));
om_b=0.09;


%% v
% Mean of \sigma_i:
X=[countryFE,yearFE(:,1:end-1)];
u=d0.y-X*(X\d0.y);
meanSigma_i=mean(sqrt(((u.^2)'*countryFE)./sum(countryFE)));

s_v=3;
d_v=2*std(avgGR)/meanSigma_i;


%% \alpha
pI1=(d0.I'*countryFE)./sum(countryFE);
pI1=max(0.05,min(pI1,0.95));
log_odds=log(pI1./(1-pI1));
a0=mean(log_odds);
pI1A=(((1-d0.D).*d0.I)'*countryFE)./((1-d0.D)'*countryFE);
pI1A=max(0.05,min(pI1A(isfinite(pI1A)),0.95));
log_oddsA=log(pI1A./(1-pI1A));
pI1D=((d0.D.*d0.I)'*countryFE)./(d0.D'*countryFE);
pI1D=max(0.05,min(pI1D(isfinite(pI1D)),0.95));
log_oddsD=log(pI1D./(1-pI1D));
yA=(((1-d0.D).*d0.y)'*countryFE)./((1-d0.D)'*countryFE);
yA=yA(isfinite(yA));
yD=((d0.D.*d0.y)'*countryFE)./(d0.D'*countryFE);
yD=yD(isfinite(yD));
om_a=sqrt(var(log_oddsA)-var(yA));


%% f
f0=0;
om_f=sqrt(max(0.01,var(log_oddsD)-var(yD)-om_a^2));


%% \varsigma
s_vs=3;
d_vs=2*om_f/4;


%% 
clearvars -except a0 b0A b0D d_v d_vs f0 om_a om_b om_f om_th s_v s_vs th0


%%
save('prior_robustTH.mat')



